% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Code for simulations based on VAR(2) specification
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; clc; close all;
rng(12345);
addpath('functions');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compile C code
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mex 'functions/buildVAR.c'
mex 'functions/buildForwards.c'
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Settings
%   3-factor model: factorFwdMats = [1,5,10];
%   5-factor model: factorFwdMats = [1,3,5,7,10];
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
o.dataSet               = 'gsw';
o.freq                  = 'm';
o.startDate             = 197201;
o.endDate               = 202212;
p.N                     = 120;
p.yldMats                 = (12:12:120);
p.factorFwdMats           = [1,5,10];
p.BH_B_BC               = 0;
p.CG_B_BC               = 0;
sp.S                    = 5000;
sp.B                    = 399;
sp.T                    = 600;
sp.M                    = 'ROT';

p.K                     = numel(p.factorFwdMats);
p.H                     = p.yldMats(1);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load bond data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' '); disp('Loading Data...'); disp(' ');
%
dataAll = loadBondData(p.N);
dataTag = ['dataAll.', o.dataSet, '.', o.freq];
for i = {'prc', 'yld', 'ret', 'xret', 'fwd', 'dret'}, eval([i{1}, 'Raw = ', dataTag, '.', i{1}, ';']);  end
data.mats = (1:p.N);
datesBonds = eval([dataTag, '.dates']);
data.period = eval([dataTag, '.period']);
datesCommon = datesBonds;
%%
if ~isempty(o.startDate), datesCommon = datesCommon(datesCommon >= o.startDate); end
if ~isempty(o.endDate), datesCommon = datesCommon(datesCommon <= o.endDate); end
%
dates = datesCommon;
idxBonds = ismember(datesBonds,dates);
%
data.prc = prcRaw(idxBonds,1:p.N);
data.yld = yldRaw(idxBonds,1:p.N);
data.fwd = fwdRaw(idxBonds,1:p.N);
data.dret = dretRaw(idxBonds,1:p.N);
data.ret = retRaw(idxBonds,1:p.N);
data.xret = xretRaw(idxBonds,1:p.N);
data.T = size(data.prc,1);
% 
data2.yld = data.yld(:,p.yldMats);
data2.fwd = [data2.yld(:,1) data2.yld(:,2:end).*(ones(data.T,1)*p.yldMats(2:end))/data.period-data2.yld(:,1:end-1).*(ones(data.T,1)*p.yldMats(1:end-1))/data.period];
data2.dret = zeros(size(data2.yld));
data2.dret(1:p.H,:) = NaN;
data2.dret(p.H+1:end,2:end) = data2.fwd(1:end-p.H,2:end) - data2.fwd(p.H+1:end,1:end-1);
data2.xret = cumsum(data2.dret,2);
%
if strcmp(sp.T,'SAME'), sp.T=data.T; end
if strcmp(sp.M,'ROT'), sp.M=ceil((sp.T*p.N)^(2/5)); end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LLSW data-driven truncation parameters
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sp.nu_LLSW_EWC = floor(0.4*(sp.T-p.H)^(2/3));
sp.S_LLSW_NW = ceil(1.3*sqrt(sp.T-p.H));

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate parameters of data generating process
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = data2.fwd(:,p.factorFwdMats);
X1 = [ones(data.T,1) X];
X12 = [ones(data.T,1) [X(2:end-1,:); NaN(2,p.K)] [X(1:end-1,:); NaN(1,p.K)]];
PsiHat = inv(X12(1:end-2,:)'*X12(1:end-2,:))*X12(1:end-2,:)'*X(3:end,:);
PhiHat = [PsiHat(2:end,:)'; eye(p.K) zeros(p.K)];
muHat =  PsiHat(1,:)';
VHat = X(3:end,:) - X12(1:end-2,:)*PsiHat;
SigHat = VHat'*VHat/(data.T - numel(PsiHat));
loadHat = inv(X1'*X1)*X1'*data2.fwd;
eHat = data2.fwd - X1*loadHat;
tmp = eHat'*eHat/(data.T-numel(loadHat));
[tmpS,tmpV,tmpD] = svd(tmp);
SigEHatSqrt = tmpS*((tmpV).^(1/2))*tmpD';
data2.SigEHatSqrt = SigEHatSqrt;
data2.SigHat = SigHat;
data2.PsiHat = PsiHat;
data2.PhiHat = PhiHat;
data2.muHat = muHat;
data2.X = X;
data2.X1 = X1;
data2.X12 = X12;
data2.loadHat = loadHat;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate in-sample regressions
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
betaHatAllData = NaN(numel(p.yldMats),size(X12,2));
for i = 1:numel(p.yldMats)
    % (1:end-2) and (3:end) is 1-step ahead so that p.H steps ahead is as given
    betaHatAllData(i,:) = regress(data2.xret(p.H+2:end,i),X12(1:end-p.H-1,:));
end

betaHatAll = NaN(sp.S,numel(p.yldMats),size(X12,2));
betaSEHatAll = NaN(sp.S,numel(p.yldMats),size(X12,2));
bsBetaHatAll_CG = NaN(sp.S,sp.B,numel(p.yldMats),size(X12,2));
bsBetaHatAll_BH = NaN(sp.S,sp.B,numel(p.yldMats),size(X12,2));
bsBetaHatAll_FullOracle = NaN(sp.S,sp.B,numel(p.yldMats),size(X12,2));
bsBetaSEHatAll_CG = NaN(sp.S,sp.B,numel(p.yldMats),size(X12,2));
bsBetaSEHatAll_BH = NaN(sp.S,sp.B,numel(p.yldMats),size(X12,2));
bsBetaSEHatAll_FullOracle = NaN(sp.S,sp.B,numel(p.yldMats),size(X12,2));
bsBetaSEHatAll_LLSW_NW = NaN(sp.S,numel(p.yldMats),size(X12,2));
bsBetaSEHatAll_LLSW_EWC = NaN(sp.S,numel(p.yldMats),size(X12,2));
simEigvals = NaN(sp.S,numel(p.yldMats));

parpool;
fprintf('\nProgress:\n');
fprintf(['\n' repmat('.',1,sp.S) '\n\n']);

tic
parfor s=1:sp.S
    simOut = simsVAR2Fcn(s,data2,data,p,sp);
    betaHatAll(s,:,:) = simOut.betaHatAll;
    betaSEHatAll(s,:,:) = simOut.betaSEHatAll;
    bsBetaHatAll_CG(s,:,:,:) = simOut.bsBetaHatAll_CG;
    bsBetaHatAll_BH(s,:,:,:) = simOut.bsBetaHatAll_BH;
    bsBetaHatAll_FullOracle(s,:,:,:) = simOut.bsBetaHatAll_FullOracle;
    bsBetaSEHatAll_CG(s,:,:,:) = simOut.bsBetaSEHatAll_CG;
    bsBetaSEHatAll_BH(s,:,:,:) = simOut.bsBetaSEHatAll_BH;
    bsBetaSEHatAll_FullOracle(s,:,:,:) = simOut.bsBetaSEHatAll_FullOracle;
    bsBetaSEHatAll_LLSW_NW(s,:,:) = simOut.betaSEHatAll_LLSW_NW;
    bsBetaSEHatAll_LLSW_EWC(s,:,:) = simOut.betaSEHatAll_LLSW_EWC;
    simEigvals(s,:) = simOut.simEigvals;
    fprintf('\b|\n');
end
disp(' ');
toc
disp(' ');
delete(gcp('nocreate'));

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% True parameters
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Et = eye(numel(p.yldMats));
Et(1,:) = [];
Eth = eye(numel(p.yldMats));
Eth(end,:) = [];

tmp = 0;
for i = 1:p.H, tmp = tmp + PhiHat^(i-1)*[muHat; zeros(p.K,1)]; end

beta0 = NaN(size(squeeze(mean(betaHatAll(:,:,1:end)))));
beta0(2:end,2:end) = cumsum(Et*[loadHat(2:end,:)' zeros(size(loadHat(2:end,:)'))] - Eth*loadHat(2:end,:)'*[eye(3) zeros(3)]*PhiHat^p.H);
beta0(2:end,1) = cumsum((Et-Eth)*loadHat(1,:)' - Eth*loadHat(2:end,:)'*[eye(3) zeros(3)]*tmp);

pval_Size_CG = NaN(sp.S,numel(p.yldMats),size(X12,2));
pval_Size_BH = NaN(sp.S,numel(p.yldMats),size(X12,2));
pval_Size_FullOracle = NaN(sp.S,numel(p.yldMats),size(X12,2));
pval_Power_CG = NaN(sp.S,numel(p.yldMats),size(X12,2));
pval_Power_BH = NaN(sp.S,numel(p.yldMats),size(X12,2));
pval_Power_FullOracle = NaN(sp.S,numel(p.yldMats),size(X12,2));
tstat_Size_NW = NaN(sp.S,numel(p.yldMats),size(X12,2));
tstat_Power_NW = NaN(sp.S,numel(p.yldMats),size(X12,2));
tstat_Size_LLSW_EWC = NaN(sp.S,numel(p.yldMats),size(X12,2));
tstat_Power_LLSW_EWC = NaN(sp.S,numel(p.yldMats),size(X12,2));
tstat_Size_LLSW_NW = NaN(sp.S,numel(p.yldMats),size(X12,2));
tstat_Power_LLSW_NW = NaN(sp.S,numel(p.yldMats),size(X12,2));

for s = 1:sp.S
    for i = 2:numel(p.yldMats)
        for j = 1:size(X12,2)
            tstat = betaHatAll(s,i,j)./betaSEHatAll(s,i,j);
            tstat0  = (betaHatAll(s,i,j)-beta0(i,j))./betaSEHatAll(s,i,j);
            tstat_Size_NW(s,i,j) = tstat0;
            tstat_Power_NW(s,i,j) = tstat;
            %
            tmp = (bsBetaHatAll_CG(s,:,i,j) - betaHatAll(s,i,j))./bsBetaSEHatAll_CG(s,:,i,j);
            pval_Size_CG(s,i,j) = mean(abs(tstat0) < abs(tmp));
            pval_Power_CG(s,i,j) = mean(abs(tstat) < abs(tmp));
            tmp = (bsBetaHatAll_BH(s,:,i,j) - betaHatAll(s,i,j))./bsBetaSEHatAll_BH(s,:,i,j);
            pval_Size_BH(s,i,j) = mean(abs(tstat0) < abs(tmp));
            pval_Power_BH(s,i,j) = mean(abs(tstat) < abs(tmp));
            tmp = (bsBetaHatAll_FullOracle(s,:,i,j) - beta0(i,j))./bsBetaSEHatAll_FullOracle(s,:,i,j);
            pval_Size_FullOracle(s,i,j) = mean(abs(tstat0) < abs(tmp));
            pval_Power_FullOracle(s,i,j) = mean(abs(tstat) < abs(tmp));                                    
            %
            tstat_Power_LLSW_EWC(s,i,j) = betaHatAll(s,i,j)./bsBetaSEHatAll_LLSW_EWC(s,i,j);
            tstat_Size_LLSW_EWC(s,i,j) = (betaHatAll(s,i,j)-beta0(i,j))./bsBetaSEHatAll_LLSW_EWC(s,i,j);
            %
            tstat_Power_LLSW_NW(s,i,j) = betaHatAll(s,i,j)./bsBetaSEHatAll_LLSW_NW(s,i,j);
            tstat_Size_LLSW_NW(s,i,j) = (betaHatAll(s,i,j)-beta0(i,j))./bsBetaSEHatAll_LLSW_NW(s,i,j);            
        end
    end
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CG and BH
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nominal Size 10%
squeeze(mean(pval_Size_CG < 0.1))
squeeze(mean(pval_Size_BH < 0.1))
squeeze(mean(pval_Size_FullOracle < 0.1))
squeeze(mean(pval_Power_CG < 0.1))
squeeze(mean(pval_Power_BH < 0.1))
squeeze(mean(pval_Power_FullOracle < 0.1))
% Nominal Size 5%
squeeze(mean(pval_Size_CG < 0.05));
squeeze(mean(pval_Size_BH < 0.05));
squeeze(mean(pval_Power_CG < 0.05));
squeeze(mean(pval_Power_BH < 0.05));

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Standard Newey-West
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
squeeze(mean(abs(tstat_Size_NW) > norminv(.95),1))
squeeze(mean(abs(tstat_Size_NW) > norminv(.975),1))

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LLSW: EWC
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nominal Size 10%
squeeze(mean(abs(tstat_Size_LLSW_EWC) > tinv(.95,sp.nu_LLSW_EWC),1))
% Nominal Size 5%
squeeze(mean(abs(tstat_Size_LLSW_EWC) > tinv(.975,sp.nu_LLSW_EWC),1))

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LLSW: NW
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 0.95 and 0.975
cv_LLSW_NW_10 = load('../input/cv_nw_10.mat');
cv_LLSW_NW_5 = load('../input/cv_nw_05.mat');
[~, idx] = min(abs(sp.S_LLSW_NW/sp.T-cv_LLSW_NW_10.cv_nw_10(:,1,1)));
% Nominal Size 10%
squeeze(mean(abs(tstat_Size_LLSW_NW) > sqrt(cv_LLSW_NW_10.cv_nw_10(idx,2,1)),1))
% Nominal Size 5%
squeeze(mean(abs(tstat_Size_LLSW_NW) > sqrt(mean(cv_LLSW_NW_5.cv_nw_05(idx,2,1))),1))
% Power Size 10%
squeeze(mean(abs(tstat_Power_LLSW_NW) > sqrt(cv_LLSW_NW_10.cv_nw_10(idx,2,1)),1))

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eigenvalue share
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mean(simEigvals(:,1)./sum(simEigvals,2))
mean(simEigvals(:,2)./sum(simEigvals,2))
mean(simEigvals(:,3)./sum(simEigvals,2))

tmpNumFac = numel(p.factorFwdMats);
if p.BH_B_BC>0, tmpBC = 'simsVAR2BC_'; else tmpBC = 'simsVAR2NoBC_'; end
save(['../output/', tmpBC, num2str(numel(p.factorFwdMats)), 'Fac_S', num2str(sp.S), '_T', num2str(sp.T)], 'pval*', 'tstat*', 'beta0', 'sp', 'simEigvals');

